library(ggplot2)
half_life <- 1.387 # half life of 10Be in Million years
initial_ratio <- 13.87e-8 # Initial Authigenic measured in Gorongosa
AMS_limit <- 1e-15 # technical limit of AMS
decay_rate <- function(t){
return(exp(-log(2) * t / half_life))
}
times <- seq(0,50,.1) # time range in Million years
decay <- decay_rate(times)
data <- data.frame(time = times, decay = decay)
x <- 19.5
y_interp <- approx(data$time, data$decay, x)$y
point_data <- data.frame(time = x, decay = y_interp)
y <- AMS_limit / initial_ratio
x_interp <- approx(data$decay, data$time, y)$y
vline_data <- data.frame(time = x_interp, decay = y)
ggplot(data, aes(x = time, y = decay)) +
geom_line(size = 2) +
geom_point(data = point_data, color = "blue", size = 4) +
geom_text(data = point_data, aes(label = round(decay,10)), vjust = -1, color = "blue") +
ggtitle("Age limit for 13.87e-8 Aut. Init. ~ 37.5 Ma") +
xlab("Time (Ma)") +
ylab("Fraction remaining") +
scale_x_continuous(breaks = seq(0,50,2)) +
scale_y_continuous(label = scales::percent) +
theme_minimal() + theme(text = element_text(size = 18)) +
geom_vline(xintercept = vline_data$time, linetype = "dashed",
colour = "red", size = 2)
ggplot(data, aes(x = time, y = decay)) +
geom_line(size = 2) +
geom_point(data = point_data, color = "blue", size = 4) +
geom_text(data = point_data, aes(label = round(decay,10)), vjust = -2, color = "blue") +
geom_text(data = point_data, aes(label = time), vjust = 2, color = "black") +
ggtitle("Plot in log, red dash age is where AMS detection limit 1e-15 intercepts") +
xlab("Time (Ma)") +
ylab("Fraction remaining (log scale)") +
scale_x_continuous(breaks = seq(0,50,2)) +
scale_y_continuous(trans = "log", label = scales::percent) +
theme_minimal() + theme(text = element_text(size = 18)) +
geom_vline(xintercept = vline_data$time, linetype = "dashed",
colour = "red", size = 2)
Authigenic plots
# Initial Authigenic values seem to range from -6 to -12 in the literature
initial <- 10^(seq(-6,-12,length.out = length(times)))
data$initial <- initial
detected <- expand.grid(initial = data$initial, decay = data$decay)
detected$detected <- detected$initial * detected$decay
# Set threshold and range for z values
threshold <- 1e-15 # AMS limit is 1 ppq
range <- 0.1e-15
# Create new data frame with x and y values that are true for z being approximately 1e-15
df <- data.frame(x = detected$initial, y = detected$decay, z = detected$detected)
df <- df[abs(df$z - threshold) <= range, c("x", "y")]
resampled_curve_indexes <- sample(1:nrow(df), 50, replace=F)
resampled_curve <- df[resampled_curve_indexes,]
ggplot(data = df, aes(x = x, y = y)) +
geom_point(colour = "red", size = 0.5) +
ggtitle("Curve of AMS detection limit (1e-15) in 10Be/9Be") +
xlab("Initial Authigenic (linear scale)") +
scale_x_continuous(breaks = c(1e-12, 1e-6)) +
ylab("Decay (linear scale)") +
scale_y_continuous(label = scales::percent) +
theme_minimal() + theme(text = element_text(size = 15))
ggplot() +
stat_summary_2d(data = detected,
aes(x = initial, y = decay, z = detected, colour = detected),
bins = 300) +
ggtitle("Remaining Isotopes in 10Be/9Be (limit = 1e-15)") +
xlab("Initial Authigenic (ratio in orders of magnitudes)")+
scale_x_continuous(trans = "log", breaks = 10^seq(-12,-6,1)) +
ylab("Decay (% in log scale)") +
scale_y_continuous(trans = "log", label = scales::percent, n.breaks = 10) +
scale_fill_viridis_c(option = "C", trans = "log", breaks = 10^seq(-15,-7,1)) +
theme_minimal() +
geom_smooth(data = resampled_curve, aes(x = x, y = y), colour = "black", method = "lm")
## `geom_smooth()` using formula 'y ~ x'
detected$time <- rep(times, each = length(decay))
limitset <- detected[detected$detected > 1e-15,]
library(plotly)
# Create a scatter plot with x, y, and z variables
p1 <- plot_ly(limitset, x = ~initial, y = ~decay, z = ~time, type = 'scatter3d',
mode = 'markers', color = ~log(detected)) %>%
layout(scene = list(xaxis = list(title = 'Initial Authigenic', type = 'linear', tickformat = ".2e"),
yaxis = list(title = 'Decay ratio (%)', type = 'linear'),
zaxis = list(title = 'Time (Million years)', type = 'linear')))
# Show the plot
p1